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Zero— temperature study of vacancies in solid 4 He 
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Abstract. We are investigating the properties of vacancies in solid He with the exact zero 
£f>l , temperature SPIGS method. Our aim is to study the possibility of phase-separation between 

vacancies and the perfect crystal. We find a significant correlation among vacancies both in 
two and in three dimensional solid 4 He. In 3D we have considered up to 5 vacancies in a hep 
crystal at a density p — 0.0293A -3 , close to the melting, and in 2D we have considered up to 6 
' vacancies in a triangular crystal at p = 0.0765A -2 . In all the considered cases the correlation 

i <~| among vacancies seems to display an exponential decay that would suggest the presence of a 

bound state. The systematically sublinear dependence on the number of vacancies found for the 
activation energy strengthens the indication of an attractive interaction, even if we have direct 
evidence that vacancies do not form a single compact cluster. 

B 

Vacancies in a quantum solid like solid Helium have attracted a great interest because of their 
unique properties pQ : extensive theoretical and experimental work on solid Helium has indicated 
that vacancies have an high probability of tunneling from site to site due to the large zero-point 
motion of the atoms. This distinguishes helium from other solids in which the vacancies are 
localized at low temperature. Despite of this large effort, our understanding of vacancies in 
solid Helium is not complete [TJ |2j. Moreover, much of the interest in vacancies in solid 4 He 
£NJ ! resulted from theoretical speculation [3] that the vacancy concentration could be finite even at 

zero temperature (zero-point vacancies) and that this could result in a Bose-Einstein condensate 
(BEC) in the solid OH]. In that case, solid 4 He would be a super solid and it could display some 
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superfluid properties such as non classical rotational inertia (NCRI) [5]. Zero-point vacancies 
have revealed to be elusive to the experimental observation [HE]- On the contrary, NCRI effects 
O ! have been recently observed in solid 4 He [7], giving a renewed interest to the supersolid state of 

matter [8]. Even for this topic, our understanding is far from being complete 0(9]. 

Interaction among vacancies in solid 4 He have been recently studied by means of an elastic 
theory [TO]. Here we systematically investigate the properties of vacancies in solid 4 He by means 
of a full microscopic technique such as the Shadow Path Integral Ground State (SPIGS) [TT] . 
The SPIGS method allows to recover "exact" expectation values on the true ground state of 
4 He systems, and it is particularly indicated for non homogeneous and defected systems. We 
have considered vacancies both in three dimensional (3D) and in two dimensional (2D) solid 
4 He. The study of vacancy properties in 3D solid 4 He at T = OK is relevant for the possible 
supersolid phase. In fact, if vacancies form a tight-bound state, as recently suggested by finite 
temperature PIMC simulations [12], they would separate breaking the mechanism that would 
induce BEC [13] . The 2D system represents a good reference model for solid 4 He adsorbed 
on planar substrates, such as graphite, which are under experimental investigation, and it 
is also relevant for its conceptual counterpart: the Abrikosov lattice of flux lines in type-II 
superconductors [H]. 

Dealing with low temperature properties, 4 He atoms are described as structureless zero-spin 



bosons, interacting through a realistic two-body potential, that we assume to be the HFDHE2 
Aziz potential |15j . The Path Integral Ground State (PIGS) [16] method recovers more and 
more accurate approximations of the true ground state wave function by successive short time 
projections in imaginary time of a trial wave function. The true ground state expectation values 
will be reached only in the limit of infinite imaginary time; in practice, being the convergence 
exponential [T7], after a suitable number of small time projections the PIGS values differ from to 
the true ones within the statistical errors and further projections affect no more the computed 
values. In this sense, the results provided by PIGS are exact. Moreover we have recently shown 
that the PIGS method is unaffected by any variational bias due to the choice of trial wave 
function [T8] . 

When a shadow wave function (SWF) [19, 20J is employed as trial wave function, the SPIGS 
method is recovered [11 . The use of SWF is justified both by the fact that it gives the best 
available variational description of solid and liquid 4 He |20| and by the fact that no equilibrium 
positions for the solid phase are requested allowing to easily describe also disordered systems |llj . 
We have employed the pair-product approximation for the imaginary time projector [21j and 
the imaginary time step 5t = 1/40K -1 has been chosen in order to ensure a good accuracy and 
a reasonable computational effort [18] . 

In order to calculate the activation energies of vacancies, we perform a number of independent 
simulations: a simulation for the perfect crystal and a simulation of the crystal with n vacancies 
for different n values. In presence of vacancies, after removing n particles, we rescale the 
dimensions of the simulation box to reset the system to the original density. This is performed 
to circumvent the need of correcting the energy due to density change caused by the inclusion 
of vacancies. Thus the activation energy of n vacancies in a crystal with M lattice sites is 
AE n = [e(M — n) — e(M)](M — n) where e(iV) is the energy per particle for the system containing 
iV particles. Starting from an initial configuration corresponding to a perfect lattice in which 
n particles have been removed, we find that, in all the considered cases, the crystalline state 
is stable and the periodicity is such that the number of maxima in the local density remains 
equal to M. This means that we have a crystal with n mobile vacancies. When dealing with the 
three dimensional solid 4 He we have considered a box with periodic boundary conditions (pbc) 
built to fit an hep crystal at p = 0.0293A~ 3 with M = 180 lattice positions. We have computed 
the activation energy for up to 5 vacancies. We find AE\ = 15.4 ± 0.8K, AE2 = 30.6 ± 0.8K, 
AE 3 = 43.9 ± 0.8, AEi = 53.7 ± 0.8K and AE 5 = 63.4 ± 0.7K. We can see that the dependence 




Figure 1. Vacancy-vacancy corre- 
lation function g vv (r) computed in 
a hep crystal with M = 720 lattice 
positions and n = 3 vacancies at the 
density p = 0.0293A" 3 for two dif- 
ferent imaginary projection time r. 
g vv {r) is normalized with the corre- 
lation function of 3 vacancies ran- 



domly distributed on the same lat- 



of AE n from n is systematically sublinear, suggesting the existence of some kind of attractive 
interaction among vacancies. A similar result is obtained also in the two dimensional case. We 
have considered a box with pbc housing a triangular crystal at p = 0.0765A -2 with M = 240 
lattice positions and a number of vacancies from 1 to 6. The obtained activation energies 
are A#i = 7.0 ± 0.3K, AE 2 = 12.2 ± 0.5K, AE 3 = 15.5 ± 0.3, AE A = 19.1 ± 0.3K and 
AEq = 25.1 ± 0.4K. Even in this case, the dependence of AE n on n is systematically sublinear. 

In order to infer whether this interaction is strong enough to give rise to a bound state and 
to a phase-separation we have computed also a vacancy-vacancy correlation function g vv (r) 
by histogramming the relative distances of the vacancies during the Monte Carlo sampling. 
The determination of the vector positions of the vacancies in a crystalline configurations is far 
from being trivial due to the large zero point motion, to high vacancy mobility and because 
in our algorithm the center of mass is not fixed. The first step of our analysis is, given a 
configuration explored by the Metropolis sampling, to find the crystal lattice {R±} which best 
fits the positions of the real particles {fj}. This is achieved by optimizing the position of the 
center of mass of the reference crystal lattice in order to maximize the Gaussian local density 

Y^jLi Y^f=i ex P[ — a (^j — ^|°^) 2 ]> where a is a parameter suitably chosen in order to have the 
best efficiency. The position of a vacancy is then recovered by coarse graining procedure, which 
allocates in subsequent steps each particle to an exclusive lattice position by considering larger 
and larger distances, up to a maximum distance 2a (a is the lattice parameter). If it turns out 
impossible to allocate all the particles within a distance of 2a from a reference lattice position, 
the configuration is discarded. In order to ensure the convergence of g vv (r) we have considered 
two different starting configurations: one with the n vacancies at random lattice positions in the 
crystal, and one with the n vacancies in a compact cluster configuration. 

Our results for 3 vacancies in a 3D hep crystal with M = 720 lattice positions at p = 
0.0293A- 3 are plotted in Fig. [T] for two different r values. The large values of g vv (r) found 
in the short distance range enforce the suggestion of a strong attractive interaction among 
vacancies, moreover g vv (r) shows an exponential decay for distance up to about 15A. Despite 
of the substantial depletion at large distances, vacancies are able to explore the whole available 
distance range, as shown by the large distance plateau of g vv {r). Our present data seem to be 
not compatible with the phase-separation picture. However, whereas up to 10A g vv (r) seems at 
convergence with respect to the r value, this is not true for the tails where g vv {r) still depends 
on r (see Fig. [TJ. Computations with larger r are underway. 
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Figure 2. Vacancy- vacancy cor- 
relation function g vv {r) computed 
in a triangular crystal with M = 
240 lattice positions and n = 2, 
3 and 4 vacancies at the den- 
sity p = 0.0765A" 2 . The total 
projection imaginary time is r = 
0.775K -1 . Correlation functions 
computed with larger r, not shown 
here, do not show significant depar- 
tures from the plotted ones. g vv (r) 
is normalized with the correlation 
function of n vacancies randomly 
distributed on the same lattice. 



Similar results are found also in the two dimensional case. Our results for 2, 3 and 4 
vacancies with a total projection time r = 0.775K -1 (larger than in the 3D case) are reported in 
Fig. [2] Contrarily to the 3D results, no large distance plateaus are found in g vv (r), suggesting a 
bound state. When n = 2, g vv {r) turns out to be well described by an exponential decay with 
correlation length A ~ 3A. For larg er Qw{j'^} shows a more complex behavior, and it seems 
possible to recognize two regimes characterized by different A values. For example, when n = 4, 
A ~ 14.3A for r < 12A and A ~ 2.97A for r > 12A. This large r behavior could be significantly 
affected by the pbc, simulations with larger boxes are underway. It should be noticed that the 
same g vv (r) is recovered both when the simulations starts with the vacancies placed in a compact 
cluster configuration and randomly placed in the crystalline lattice. Thus, even if initially placed 
in a compact cluster configuration vacancies explore the whole available distance range: this is 
a direct evidence that vacancies do not form a compact cluster. As already pointed out, our 
algorithm discards some configurations where it fails to allocate the particles to the reference 
lattice sites within the maximum distance of 2a. By looking at these configurations in the 2D 
case for the larger n considered values, we can see that sometime the system rearranges the 
point defects into couples of dislocations, which in 2D are still point defects. The appearance of 
such topological excitations is observed also in classical 2D solids [22] . 

To conclude, we have presented a study of vacancies in solid 4 He in 3D and 2D crystals. 
We have computed the activation energy AE n for different numbers of vacancies n and we 
have found a systematic sublinear dependence of AE n on n, which gives the indication of an 
attractive interaction between point defects. This is quite similar to what is found for classical 
solids [23]. Following the positions of the vacancies during the Monte Carlo sampling we were 
able to compute a vacancy- vacancy correlation function g vv (r). For all the considered cases 
g vv (r) seems to display an exponential decay which suggests that vacancies form a bound-state. 
But, at the same time, vacancies are always able to explore the whole available distance range 
hence do not form a compact cluster. Further investigations are needed in order to have more 
quantitative details on vacancy-vacancy interactions. 
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